!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

   SUBROUTINE TRIANGLE_GRID_EDGE

! NEEDS UPDATE TO FIX ERROR MESSAGES

!==============================================================================!
!  This program is used to define the non-overlapped, unstructured             !
!  triangular meshes used for flux computations. The mesh could be             !
!  created using the commerical software called "sms8.0" or other              !
!  mesh generation programs. The mesh file generated by sms8.0 can             !
!  be directly used for this subroutine, while the mesh file                   !
!  generated using other programs must be coverted the data format             !
!  to meet the required format used here.                                      !
!==============================================================================!
!     variable list:							       !
!  vx(m)    :: vx(i) = x-coordinate of node i (input from mesh)	               !
!  vy(m)    :: vy(i) = y-coordinate of node i (input from mesh)	               !
!  nv(n,3)  :: nv(i:1-3) = 3 node numbers of element i                         !
!  xc(n)    :: xc(i) = x-coordinate of element i (calculated from vx)          !
!  yc(n)    :: yc(i) = y-coordinate of element i (calculated from vy)          !
!  cor(n)   :: cor(i) = f plane coriolis at element i                          !
!                                                                              !
!  nbe(n,3) :: nbe(i,1->3) = element index of 1->3 neighbors of element i      !
!  isbce(n) :: flag if element is on the boundary, see below for values        !
!  isonb(m) :: flag is node is on the boundary, see below for values           !
!                                                                              !
!  ntve(m)  :: the number of neighboring elements of node m                    !
!  nbve(m,ntve(m)) :: nbve(i,1->ntve(i)) = ntve elements containing node i     !
!  nbvt(m,ntve(m)) :: nbvt(i,j) = the node number of node i in element         ! 
!                     nbve(i,j) (has a value of 1,2,or 3)                      ! 
!                                                                              !
!   ne       ::   number of unique element edges                               !
!  iec(ne,2) :: nbe(i,1->2) cell number of cell(s) connected to edge i         !
!   isbc(ne) :: flag marking edge property                                     !
!              isbc(i) = 0:  element edge i is in the interior                 !
!              isbc(i) = 1:  element edge i is on the boundary                 !
!ienode(ne,2):: ienode(i,1->2) node numbers at each end of element edge i      ! 
!  xijc(ne)  :: xijc(i) = x-coordinate of mid point of element edge i          ! 
!  yijc(ne)  :: yijc(i) = y-coordinate of mid point of element edge i          ! 
!  dltxyc(ne):: dltxyc(i) = length of element edge i                           !
!  dltxc(ne) :: dltxc(i) = deltax (x-projection) of element edge i             !
!  dltyc(ne) :: dltyc(i) = deltay (y-projection) of element edge i             !
!  sitac(ne) :: sitac(i) =  arctg(dltyc,dltxc) (angle of inclination of edge)  !
!                                                                              !
!==============================================================================!
!     classification of the triangles nodes, and edges                         !
!                                                                              !
!     isonb(i)=0:  node in the interior computational domain                   !
!     isonb(i)=1:  node on the solid boundary                                  !
!     isonb(i)=2:  node on the open boundary                                   !
!                                                                              !
!     isbce(i)=0:  element in the interior computational domain                !
!     isbce(i)=1:  element on the solid boundary                               !
!     isbce(i)=2:  element on the open boundary                                !
!     isbce(i)=3:  element with 2 solid boundary edges                         !
!                                                                              !
!      isbc(i)=0:  element edge in interior                                    !
!      isbc(i)=1:  element edge on boundary                                    !
!==============================================================================!


!==============================================================================|
!   FIND NEIGHBORING ELEMENTS, MARK BOUNDARY NODES AND ELEMENTS                |
!									       |
!   NBE(N,3) :  NBE(I,J) = IDENTITY OF NBOR ELMNT TO TRI I ON EDGE J           |
!   IBCE(N)  :  DESCRIBED IN SUBROUTINE HEADING			               |	
!   ISONB(M):  DESCRIBED IN SUBROUTINE HEADING			               |	
!==============================================================================|
   USE ALL_VARS
   USE MOD_SPHERICAL


   USE MOD_PAR

   USE MOD_OBCS
   IMPLICIT NONE
   INTEGER, ALLOCATABLE, DIMENSION(:,:) :: TEMP,TEMP2,NB_TMP,ISET
   INTEGER I,J,DIF1,DIF2,DIF3,II,JJ,NTMP,NCNT,INEY,NFLAG
   INTEGER ITMP1,ITMP2,ITMP3,JN,JJB,IBCETMP,NCTMP,NCETMP,NPT
   INTEGER, ALLOCATABLE :: CELLS(:,:),CELLCNT(:),NBET(:,:)
   REAL(SP)  DTMP
   INTEGER N1,N2,N3,J1,J2,J3,SBUF, IERR
   
   ! STORAGE FOR ISONB, ISBCE TO RUN EXCHANGE AND SET HALO'S
   REAL(SP), ALLOCATABLE :: FTEMP(:)

   Character(len=8)  :: tmpstr
   Character(len=200):: errstr


# if defined (SPHERICAL)
   REAL(DP) XXC,YYC,side,&
            DELTX,DELTY,x1_dp,y1_dp,x2_dp,y2_dp,DTMP_dp
   REAL(DP) XTMP,XTMP1	    
# else
   REAL(SP) DELTX,DELTY,ALPHA1,ALPHA2
# endif

!----------------------------REPORT--------------------------------------------!
   IF (DBG_SET(DBG_LOG))THEN
      WRITE(IPT,*  )'!'
      WRITE(IPT,*)'!           SETTING UP TRIS/ELEMENTS/EDGES/CVS          '
      WRITE(IPT,*  )'!'
   END IF
!----------------------------INITIALIZE----------------------------------------!
   
   ISBCE = 0
   ISONB = 0
   NBE   = 0

!
!----DETERMINE NBE(i=1:n,j=1:3): INDEX OF 1 to 3 NEIGHBORING ELEMENTS----------!
!
   ALLOCATE(NBET(NT,3)) ; NBET = 0
   ALLOCATE(CELLS(MT,50)) ; CELLS = 0
   ALLOCATE(CELLCNT(MT))  ; CELLCNT = 0
   DO I=1,NT
     N1 = NV(I,1) ; CELLCNT(N1) = CELLCNT(N1)+1
     N2 = NV(I,2) ; CELLCNT(N2) = CELLCNT(N2)+1
     N3 = NV(I,3) ; CELLCNT(N3) = CELLCNT(N3)+1
     CELLS(NV(I,1),CELLCNT(N1)) = I
     CELLS(NV(I,2),CELLCNT(N2)) = I
     CELLS(NV(I,3),CELLCNT(N3)) = I
   END DO
   if(maxval(cellcnt) > 50)write(ipt,*)'bad',maxval(cellcnt)
   DO I=1,NT
     N1 = NV(I,1)
     N2 = NV(I,2)
     N3 = NV(I,3)
     DO J1 = 1,CELLCNT(N1) 
     DO J2 = 1,CELLCNT(N2) 
       IF((CELLS(N1,J1) == CELLS(N2,J2)).AND. CELLS(N1,J1) /= I)NBE(I,3) = CELLS(N1,J1)
     END DO
     END DO
     DO J2 = 1,CELLCNT(N2) 
     DO J3 = 1,CELLCNT(N3) 
       IF((CELLS(N2,J2) == CELLS(N3,J3)).AND. CELLS(N2,J2) /= I)NBE(I,1) = CELLS(N2,J2)
     END DO
     END DO
     DO J1 = 1,CELLCNT(N1) 
     DO J3 = 1,CELLCNT(N3) 
       IF((CELLS(N1,J1) == CELLS(N3,J3)).AND. CELLS(N1,J1) /= I)NBE(I,2) = CELLS(N3,J3)
     END DO
     END DO
   END DO
   DEALLOCATE(CELLS,CELLCNT)

   ! OLD - USED TO MAKE GLOBAL INDEX ARRAYS TO SAVE PARALLEL OUTPUT
   ! REPLACED IN MOD_NCDIO WITH MORE GENERIC TYPE GRID DATA!
!   IF (PAR) THEN
!      ALLOCATE(NBEGL(0:NT,3)); NBEGL      = 0
!      DO I = 1,NT
!         NBEGL(I,:) = EGID_X(NBE(I,:))
!      END DO
!   ELSE
!      NBEGL => NBE
!   END IF

   IF (DBG_SET(DBG_LOG))WRITE(IPT,*)  '!  NEIGHBOR FINDING      :    COMPLETE'
!
!--ENSURE ALL ELEMENTS HAVE AT LEAST ONE NEIGHBOR------------------------------!
!
   NFLAG = 0
   DO I=1,NT
     IF(SUM(NBE(I,1:3))==0)THEN 
       NFLAG = 1
       WRITE(IPT,*)'ELEMENT ',I,' AT ',XC(I),YC(I),' HAS NO NEIGHBORS'
       CALL PSTOP
     END IF
   END DO
   IF(NFLAG == 1) CALL PSTOP
     
!
!----IF ELEMENT ON BOUNDARY SET ISBCE(I)=1 AND ISONB(J)=1 FOR BOUNDARY NODES J-!
!

   DO I=1,NT 
     IF(MIN(NBE(I,1),NBE(I,2),NBE(I,3))==0)THEN    !!ELEMENT ON BOUNDARY
       ISBCE(I) = 1
       IF(NBE(I,1) == 0)THEN 
         ISONB(NV(I,2)) = 1 ; ISONB(NV(I,3)) = 1
       END IF
       IF(NBE(I,2) ==0) THEN
         ISONB(NV(I,1)) = 1 ; ISONB(NV(I,3)) = 1
       END IF
       IF(NBE(I,3) ==0) THEN
         ISONB(NV(I,1)) = 1 ; ISONB(NV(I,2)) = 1
       END IF
     END IF
   END DO
   IF (DBG_SET(DBG_LOG)) WRITE(IPT,*)  '!  ISONB SETTING         :    COMPLETE'
        
!==============================================================================|
!             DEFINE NTVE, NBVE, NBVT                                          !
!                                                                              !
! ntve(1:m):           total number of the surrounding triangles               !
!                      connected to the given node                             !
! nbve(1:m, 1:ntve+1): the identification number of surrounding                !
!                      triangles with a common node (counted clockwise)        !
! nbvt(1:m,ntve(1:m)): the idenfication number of a given node over            !
!                      each individual surrounding triangle(counted            !
!                      clockwise)                                              !
! ntsn(1:m):           total number of surrounding nodes                       !
! nbsn(1:m, ntsn):     the identification number of surrounding nodes          !
!                      (counted clockwise)                                     !
! nbse(1:m,2*ntsn):    the identification number of control volume s           !
!                      edges between two neighbor nodes                        !
!==============================================================================|

!
!----DETERMINE MAX NUMBER OF SURROUNDING ELEMENTS------------------------------!
!
   MX_NBR_ELEM = 0
   DO I=1,M
     NCNT = 0
     DO J=1,NT
!       IF( (NV(J,1)-I)*(NV(J,2)-I)*(NV(J,3)-I)==0) NCNT = NCNT + 1
!       IF( (NV(J,1)-I) ==0 .OR. (NV(J,2)-I) ==0 .OR. (NV(J,3)-I)==0) NCNT = NCNT + 1
       IF( FLOAT(NV(J,1)-I)*FLOAT(NV(J,2)-I)*FLOAT(NV(J,3)-I) == 0.0_SP) &
         NCNT = NCNT + 1
     END DO
     MX_NBR_ELEM = MAX(MX_NBR_ELEM,NCNT)
   END DO
!   WRITE(IPT,*) 'MAXIMUM NUMBER OF NEIGHBOR ELEMENTS',MX_NBR_ELEM

# if defined(MULTIPROCESSOR)
   IF(PAR) THEN
      SBUF = MX_NBR_ELEM
      CALL MPI_ALLREDUCE(SBUF,MX_NBR_ELEM,1,MPI_INTEGER,MPI_MAX,MPI_FVCOM_GROUP,IERR)
   END IF
# endif

!
!----ALLOCATE ARRAYS BASED ON MX_NBR_ELEM--------------------------------------!
! 
   ALLOCATE(NBVE(M,MX_NBR_ELEM+1)); NBVE = 0
   ALLOCATE(NBVT(M,MX_NBR_ELEM+1)); NBVT = 0
!   ALLOCATE(NBSN(M,MX_NBR_ELEM+2))
   ALLOCATE(NBSN(M,MX_NBR_ELEM+3)); NBSN = 0
!
!--DETERMINE NUMBER OF SURROUNDING ELEMENTS FOR NODE I = NTVE(I)---------------!
!--DETERMINE NBVE - INDICES OF NEIGHBORING ELEMENTS OF NODE I------------------!
!--DETERMINE NBVT - INDEX (1,2, or 3) OF NODE I IN NEIGHBORING ELEMENT---------!
!
       
   DO I=1,M
     NCNT=0
     DO J=1,NT
!       IF( (NV(J,1)-I) == 0 .OR.  (NV(J,2)-I) == 0 .OR. (NV(J,3)-I) == 0)THEN 
        IF (FLOAT(NV(J,1)-I)*FLOAT(NV(J,2)-I)*FLOAT(NV(J,3)-I) == 0.0_SP)THEN
         NCNT = NCNT+1
         NBVE(I,NCNT)=J
         IF((NV(J,1)-I) == 0) NBVT(I,NCNT)=1
         IF((NV(J,2)-I) == 0) NBVT(I,NCNT)=2
         IF((NV(J,3)-I) == 0) NBVT(I,NCNT)=3
       END IF
     ENDDO
     NTVE(I)=NCNT
   ENDDO
!
!--Reorder Order Elements Surrounding a Node to Go in a Cyclical Procession----!
!--Determine NTSN  = Number of Nodes Surrounding a Node (+1)-------------------!
!--Determine NBSN  = Node Numbers of Nodes Surrounding a Node------------------!
!
!Lettmann&JQI   ALLOCATE(NB_TMP(M,MX_NBR_ELEM+1))
   ALLOCATE(NB_TMP(MX_NBR_ELEM+1,2))
   DO I=1,M
     IF(ISONB(I) == 0) THEN
       NB_TMP(1,1)=NBVE(I,1)
       NB_TMP(1,2)=NBVT(I,1)
       DO J=2,NTVE(I)+1
         II=NB_TMP(J-1,1)
         JJ=NB_TMP(J-1,2)
         NB_TMP(J,1)=NBE(II,JJ+1-INT((JJ+1)/4)*3)
         JJ=NB_TMP(J,1)
         IF((NV(JJ,1)-I) == 0) NB_TMP(J,2)=1
         IF((NV(JJ,2)-I) == 0) NB_TMP(J,2)=2
         IF((NV(JJ,3)-I) == 0) NB_TMP(J,2)=3
       ENDDO

       DO J=2,NTVE(I)+1
         NBVE(I,J)=NB_TMP(J,1)
       ENDDO

       DO J=2,NTVE(I)+1
         NBVT(I,J)=NB_TMP(J,2)
       ENDDO

       NTMP=NTVE(I)+1
       IF(NBVE(I,1) /= NBVE(I,NTMP)) THEN
          PRINT*, MYID,ngid(I),nTMP,'NBVE(I,nTMP) NOT CORRECT!!'
          PRINT*, "NBVE(I,:)=",EGID(NBVE(I,:))
          CALL PSTOP
       ENDIF

       IF(NBVT(I,1) /= NBVT(I,NTMP)) THEN
          PRINT*, ngid(I),'NBVT(I) NOT CORRECT!!'
          PRINT*, "NBVT(I,:)=",NBVT(I,:)
          CALL PSTOP
       END IF

       NTSN(I)=NTVE(I)

       DO J=1,NTSN(I)
         II=NBVE(I,J)
         JJ=NBVT(I,J)
         NBSN(I,J)=NV(II,JJ+1-INT((JJ+1)/4)*3)
       ENDDO

       NTSN(I)=NTSN(I)+1
       NBSN(I,NTSN(I))=NBSN(I,1)

     ELSE 
           JJB=0

       DO J=1,NTVE(I)
         JJ=NBVT(I,J)
         IF(NBE(NBVE(I,J),JJ+2-INT((JJ+2)/4)*3) == 0) THEN
           JJB=JJB+1
           NB_TMP(JJB,1)=NBVE(I,J)
           NB_TMP(JJB,2)=NBVT(I,J)
         END IF
       ENDDO

       IF(JJB /= 1) THEN
         WRITE(IPT,*) 'ERROR IN ISONB !,I,J', I,J
         CALL FATAL_ERROR("ERROR IN TGE DETERMINING ISONB")
       END IF

       DO J=2,NTVE(I)
         II=NB_TMP(J-1,1)
         JJ=NB_TMP(J-1,2)
         NB_TMP(J,1)=NBE(II,JJ+1-INT((JJ+1)/4)*3)
         JJ=NB_TMP(J,1)
         IF((NV(JJ,1)-I) == 0) NB_TMP(J,2)=1
         IF((NV(JJ,2)-I) == 0) NB_TMP(J,2)=2
         IF((NV(JJ,3)-I) == 0) NB_TMP(J,2)=3
       ENDDO

       DO J=1,NTVE(I)
         NBVE(I,J)=NB_TMP(J,1)
         NBVT(I,J)=NB_TMP(J,2)
       ENDDO

       NBVE(I,NTVE(I)+1)=0
       NTSN(I)=NTVE(I)+1
       NBSN(I,1)=I

       DO J=1,NTSN(I)-1
         II=NBVE(I,J)
         JJ=NBVT(I,J)
         NBSN(I,J+1)=NV(II,JJ+1-INT((JJ+1)/4)*3)
       ENDDO

       J=NTSN(I)
       II=NBVE(I,J-1)
       JJ=NBVT(I,J-1)
       NBSN(I,J+1)=NV(II,JJ+2-INT((JJ+2)/4)*3)
       NTSN(I)=NTSN(I)+2
       NBSN(I,NTSN(I))=I
     END IF
   END DO
   DEALLOCATE(NB_TMP)
   IF(MX_NBR_ELEM+3 -MAXVAL(NTSN) < 0)THEN
      WRITE(IPT,*)'CHECK NTSN/NBSN',MAXVAL(NTSN),MX_NBR_ELEM+3
      CALL PSTOP
   END IF
  

   ! OLD - USED TO MAKE GLOBAL INDEX ARRAYS TO SAVE PARALLEL OUTPUT
   ! REPLACED IN MOD_NCDIO WITH MORE GENERIC TYPE GRID DATA!
!   IF(PAR) THEN
!      ALLOCATE(NBSNGL(M,MX_NBR_ELEM+3))
!      ALLOCATE(NBVEGL(M,MX_NBR_ELEM+1))
!      DO I = 1,M
!         NBVEGL(I,:)=EGID_X(NBVE(I,:))
!
!         NBSNGL(I,:) = NGID_X(NBSN(I,:))
!      END DO
!   ELSE
!      NBVEGL => NBVE
!
!      NBSNGL =>NBSN
!   END IF

   IF (DBG_SET(DBG_LOG))WRITE(IPT,*)  '!  NBVE/NBVT             :    COMPLETE'
     

!==============================================================================!
!  Define the parameters of each triangular edge                               !
!                                                                              !
!  ne           :    number of unique element edges                            !
!  iec(1:ne,1:2):    counting number identifying two connected cells           !
!  isbc(1:ne):       0: triangle s edge in the interior                        !
!                    1: triangle s edge on the boundary                        !
!  ienode(1:ne,1:2): the identification number of two end points of a          !
!                    edge                                                      !
!  xijc(1:ne):       the x-coordinate location of the middle points            !
!                    of a edge                                                 !
!  yijc(1:ne):       the y-coordinate location of the middle points            !
!                    of a edge                                                 !
!  dltxyc(1:ne):     length of the edge                                        !
!  dltxc(1:ne):      vx(ienode(i,2))-vx(idnode(i,1))                           !
!  dltyc(1:ne):      vy(ienode(i,2))-vy(idnode(i,1))                           !
!  sitac(1:ne):      arctg(dltyc,dltxc)                                        !
!==============================================================================!

   ALLOCATE(ISET(NT,3),TEMP((NT)*3,2),TEMP2((NT)*3,2))
   ISET = 0
   NE = 0
   TEMP = 0 
   TEMP2 = 0
   DO I=1,NT
     DO J=1,3
       IF(ISET(I,J) == 0)THEN
         NE   = NE + 1
         INEY = NBE(I,J)
         ISET(I,J) = 1
         DO JN=1,3
           IF(I == NBE(INEY,JN)) ISET(INEY,JN) = 1
         END DO
         TEMP(NE,1) = I ; TEMP(NE,2) = INEY
         TEMP2(NE,1) = NV(I,J+1-INT((J+1)/4)*3)
         TEMP2(NE,2) = NV(I,J+2-INT((J+2)/4)*3)
       END IF
     END DO
   END DO
   DEALLOCATE(ISET)
!
!--ALLOCATE ARRAYS REQUIRING NUMBER OF EDGES-----------------------------------!
!
   ALLOCATE(IEC(NE,2))
   ALLOCATE(IENODE(NE,2))
   ALLOCATE(XIJC(NE))
   ALLOCATE(YIJC(NE))
   ALLOCATE(DLTXYC(NE))
   ALLOCATE(DLTXC(NE))
   ALLOCATE(DLTYC(NE))
   ALLOCATE(SITAC(NE))
   ALLOCATE(ISBC(NE))


   IEC(:,1) = TEMP(1:NE,1)
   IEC(:,2) = TEMP(1:NE,2)
   IENODE(:,1) = TEMP2(1:NE,1)
   IENODE(:,2) = TEMP2(1:NE,2)


   DEALLOCATE(TEMP,TEMP2)
         
         
!
!------MARK ELEMENT EDGES THAT ARE ON THE BOUNDARY-----------------------------!
!
   ISBC = 0
   DO I=1,NE
     IF((IEC(I,1) == 0) .OR. (IEC(I,2) == 0)) ISBC(I) = 1 
   END DO

!
!------CALCULATE ELEMENT EDGE METRICS------------------------------------------!
!
#  if defined (SPHERICAL)
   DO I=1,NE
     X1_DP=VX(IENODE(I,1))
     Y1_DP=VY(IENODE(I,1))
     X2_DP=VX(IENODE(I,2))
     Y2_DP=VY(IENODE(I,2))
     CALL ARCC(X2_DP,Y2_DP,&
               X1_DP,Y1_DP,XXC,YYC)
     CALL ARCX(X1_DP,Y1_DP,X2_DP,Y2_DP,SIDE)
!     DLTXC(I) = SIDE
     XTMP  = VX(IENODE(I,2))*TPI-VX(IENODE(I,1))*TPI
     XTMP1 = VX(IENODE(I,2))-VX(IENODE(I,1))
     IF(XTMP1 > 180.0_SP)THEN
       XTMP = -360.0_SP*TPI+XTMP
     ELSE IF(XTMP1 < -180.0_SP)THEN
       XTMP = 360.0_SP*TPI+XTMP
     END IF		
!     DLTXC(I) =XTMP*TPI*COS(DEG2RAD*(VY(IENODE(I,2))+VY(IENODE(I,1)))*0.5)
     DLTXC(I) =XTMP*COS(DEG2RAD*(VY(IENODE(I,2))+VY(IENODE(I,1)))*0.5_SP)
     DLTYC(I) =(VY(IENODE(I,2))-VY(IENODE(I,1)))*TPI
     CALL ARC(X1_DP,Y1_DP,X2_DP,Y2_DP,SIDE)
     DLTXYC(I) = SIDE
     SITAC(I)=ATAN2(DLTYC(I),DLTXC(I))      !be careful
     XIJC(I)=XXC                            !
     YIJC(I)=YYC                            !
   END DO
#  else
   DO I=1,NE
     DLTXC(I) =  VX(IENODE(I,2))-VX(IENODE(I,1))
     DLTYC(I) =  VY(IENODE(I,2))-VY(IENODE(I,1))
     XIJC(I)  = (VX(IENODE(I,1))+VX(IENODE(I,2)))/2.0_SP
     YIJC(I)  = (VY(IENODE(I,1))+VY(IENODE(I,2)))/2.0_SP
     DLTXYC(I)= SQRT(DLTXC(I)**2+DLTYC(I)**2)
     SITAC(I) = ATAN2(DLTYC(I),DLTXC(I))
   END DO
#  endif
    IF (DBG_SET(DBG_LOG))WRITE(IPT,*)  '!  EDGE SETUP            :    COMPLETE'


!==============================================================================!
!  read triangular mesh parameters on open boundary :                          !
!  iobce:   number of open boundary cells.                                     !
!  isbcn:   number of open boundary nodes.                                     !
!  i_obc_e: counter number of open boundary cells                              !
!  i_obc_n: counter number of open boundary nodes                              !
!==============================================================================!

!
!----TRAVERSE  BOUNDARY NODE NUMBERS AND SET ISONB(NODE)=2---------------------!
!
   DO I=1,IOBCN
     ISONB(I_OBC_N(I))=2
   ENDDO

!
!---- SET HALO VALUES IF PAR
! 

# if defined (MULTIPROCESSOR)
   IF(PAR) THEN
      CALL AEXCHANGE(NC,MYID,NPROCS,ISONB)
   END IF
# endif

!
!----DETERMINE IF ELEMENT IS ON OPEN BOUNDARY (CONTAINS EDGE ON OPEN BOUNDARY)-!
!
   IBCETMP=0
   DO I=1,N
     ITMP1=ISONB(NV(I,1))
     ITMP2=ISONB(NV(I,2))
     ITMP3=ISONB(NV(I,3))

     IF(SUM(ISONB(NV(I,1:3))) == 4) THEN
       ISBCE(I)=2
       IBCETMP =IBCETMP+1
     ELSE IF(SUM(ISONB(NV(I,1:3))) > 4) THEN
#  if defined (MULTIPROCESSOR)
       PRINT*,'SORRY, THE BOUNDARY CELL',EGID(I),'IS NOT GOOD FOR MODEL.'
#  else
       PRINT*,'SORRY, THE BOUNDARY CELL',I,'IS NOT GOOD FOR MODEL.'
#  endif       
       PRINT*,'IT HAS EITHER TWO SIDES OF OPEN BOUNDARY OR ONE OPEN BOUNDARY'
       PRINT*,'AND ONE SOLID BOUNDARY. PLEASE CHECK AND MODIFIED IT.'
       PRINT*,'THIS MESSAGE IS IN SUBROUTINE TRIANGLE_GRID_EDGE (TGE.F)'
       PRINT*,'STOP RUNNING...'
#  if !defined (PLBC)
       CALL PSTOP
#  endif

     END IF
   END DO

   DO I=1,NT
     IF((NBE(I,1)+NBE(I,2)+NBE(I,3) == 0).AND.(ISBCE(I) /= 2)) ISBCE(I)=3
     IF((NBE(I,1)+NBE(I,2) == 0).AND.(ISBCE(I) /= 2)) ISBCE(I)=3
     IF((NBE(I,2)+NBE(I,3) == 0).AND.(ISBCE(I) /= 2)) ISBCE(I)=3
     IF((NBE(I,1)+NBE(I,3) == 0).AND.(ISBCE(I) /= 2)) ISBCE(I)=3
   ENDDO

! SET HALO VALUES IF PAR
# if defined (MULTIPROCESSOR)
   IF(PAR) THEN
      CALL AEXCHANGE(EC,MYID,NPROCS,ISBCE)
   END IF
# endif


!==============================================================================!
!  xije(1:nc,1:2):  the x coordinate locations of starting and ending          !
!                   points of the control volumes edges                        !
!  yije(1:nc,1:2):  the y coordinate locations of starting and ending          !
!                   points of the control volumes edges                        !
!  niec(1:nc,1:2):  the counting number of left and right nodes                !
!                   conected to this control volumes edge from                 !
!                   starting point to ending point                             !
!  dltxe(1:nc):     the x distance of individual edges                         !
!  dltye(1:nc)      the y distance of individual edges                         !
!  dltxye(1:nc):    the length of individual edges                             !
!  ntrg(1:nc)  :    element associated with this control volume edge           !
!==============================================================================!
   NCTMP  = 0
   NCETMP = 0 

   DO I=1,NE
     IF(ISBC(I) == 0) THEN
       IF(IEC(I,1) <= N)THEN
         NCTMP=NCTMP+1
         NPT  =NCTMP
       ELSE
         NCETMP = NCETMP + 1
         NPT    = NCETMP+(3*N)
       END IF
       XIJE(NPT,1) = XC(IEC(I,1))
       YIJE(NPT,1) = YC(IEC(I,1))
       XIJE(NPT,2) = XIJC(I)
       YIJE(NPT,2) = YIJC(I)
       NIEC(NPT,1) = IENODE(I,1)
       NIEC(NPT,2) = IENODE(I,2)
       NTRG(NPT)   = IEC(I,1)
       DLTXE(NPT)  = XIJE(NPT,2)-XIJE(NPT,1)
       DLTYE(NPT)  = YIJE(NPT,2)-YIJE(NPT,1)
#      if defined (SPHERICAL)
       X1_DP=XIJE(NPT,1)
       Y1_DP=YIJE(NPT,1)
       X2_DP=XIJE(NPT,2)
       Y2_DP=YIJE(NPT,2)
       CALL ARCC(X2_DP,Y2_DP,X1_DP,Y1_DP,XXC,YYC)
! for global case when VX change from 360 to 0 or from 0 to 360 (degree)       
       XTMP  = XIJE(NPT,2)*TPI-XIJE(NPT,1)*TPI
       XTMP1 = XIJE(NPT,2)-XIJE(NPT,1)
       IF(XTMP1 >  180.0_SP)THEN
         XTMP = -360.0_SP*TPI+XTMP
       ELSE IF(XTMP1 < -180.0_SP)THEN
         XTMP =  360.0_SP*TPI+XTMP
       END IF	 
!       DLTXE(NPT)=TPI*DLTXE(NPT)*COS(DEG2RAD*YYC)
       DLTXE(NPT)=XTMP*COS(DEG2RAD*Y1_DP)
       DLTYE(NPT)=TPI*DLTYE(NPT)
       CALL ARC(X1_DP,Y1_DP,X2_DP,Y2_DP,DTMP_DP)
       DLTXYE(NPT)=DTMP_DP
       SITAE(NPT)=ATAN2(DLTYE(NPT),DLTXE(NPT))
       DLTXNE(I,1) = DLTXE(NPT)
       DLTYNE(I,1) = DLTYE(NPT)
#      else
       DTMP        = DLTXE(NPT)*DLTXE(NPT)+DLTYE(NPT)*DLTYE(NPT)
       DLTXYE(NPT) = SQRT(DTMP)
       SITAE(NPT)  = ATAN2(DLTYE(NPT),DLTXE(NPT))
#      endif

       IF(IEC(I,2) <= N)THEN
         NCTMP=NCTMP+1
         NPT  =NCTMP
       ELSE
         NCETMP = NCETMP + 1
         NPT    = NCETMP+(3*N)
       END IF
       XIJE(NPT,1)=XC(IEC(I,2))
       YIJE(NPT,1)=YC(IEC(I,2))
       XIJE(NPT,2)=XIJC(I)
       YIJE(NPT,2)=YIJC(I)
       NIEC(NPT,1)=IENODE(I,2)
       NIEC(NPT,2)=IENODE(I,1)
       NTRG(NPT)=IEC(I,2)
       DLTXE(NPT)=XIJE(NPT,2)-XIJE(NPT,1)
       DLTYE(NPT)=YIJE(NPT,2)-YIJE(NPT,1)
#      if defined (SPHERICAL)
       X1_DP=XIJE(NPT,1)
       Y1_DP=YIJE(NPT,1)
       X2_DP=XIJE(NPT,2)
       Y2_DP=YIJE(NPT,2)
       CALL ARCC(X2_DP,Y2_DP,X1_DP,Y1_DP,XXC,YYC)
! for global case when VX change from 360 to 0 or from 0 to 360 (degree)       
       XTMP  = XIJE(NPT,2)*TPI-XIJE(NPT,1)*TPI
       XTMP1 = XIJE(NPT,2)-XIJE(NPT,1)
       IF(XTMP1 >  180.0_SP)THEN
         XTMP = -360.0_SP*TPI+XTMP
       ELSE IF(XTMP1 < -180.0_SP)THEN
         XTMP =  360.0_SP*TPI+XTMP
       END IF	 
!       DLTXE(NPT)=TPI*DLTXE(NPT)*COS(DEG2RAD*YYC)
       DLTXE(NPT)=XTMP*COS(DEG2RAD*Y1_dp)
       DLTYE(NPT)=TPI*DLTYE(NPT)
       CALL ARC(X1_DP, Y1_DP, X2_DP, Y2_DP, DTMP_DP)
       DLTXYE(NPT)=DTMP_DP
       SITAE(NPT)=ATAN2(DLTYE(NPT),DLTXE(NPT))
       DLTXNE(I,2) = DLTXE(NPT)
       DLTYNE(I,2) = DLTYE(NPT)
#      else
       DTMP=DLTXE(NPT)*DLTXE(NPT)+DLTYE(NPT)*DLTYE(NPT)
       DLTXYE(NPT)=SQRT(DTMP)
       SITAE(NPT)=ATAN2(DLTYE(NPT),DLTXE(NPT))
#      endif
     ELSE IF(ISBC(I) == 1) THEN
       IF(IEC(I,1) <= N)THEN
         NCTMP=NCTMP+1
         NPT  =NCTMP
       ELSE
         NCETMP = NCETMP + 1
         NPT    = NCETMP+(3*N)
       END IF
       IF(IEC(I,1) == 0) THEN
         PRINT*, I,'IEC(I,1)===0'
         CALL PSTOP
       END IF
       XIJE(NPT,1)=XC(IEC(I,1))
       YIJE(NPT,1)=YC(IEC(I,1))
       XIJE(NPT,2)=XIJC(I)
       YIJE(NPT,2)=YIJC(I)
       NIEC(NPT,1)=IENODE(I,1)
       NIEC(NPT,2)=IENODE(I,2)
       NTRG(NPT)=IEC(I,1)
       DLTXE(NPT)=XIJE(NPT,2)-XIJE(NPT,1)
       DLTYE(NPT)=YIJE(NPT,2)-YIJE(NPT,1)
#      if defined (SPHERICAL)
       X1_DP= XIJE(NPT,1)
       Y1_DP= YIJE(NPT,1)
       X2_DP= XIJE(NPT,2)
       Y2_DP= YIJE(NPT,2)
       CALL ARCC(X2_DP,Y2_DP,X1_DP,Y1_DP,XXC,YYC)
! for global case when VX change from 360 to 0 or from 0 to 360 (degree)       
       XTMP  = XIJE(NPT,2)*TPI-XIJE(NPT,1)*TPI
       XTMP1 = XIJE(NPT,2)-XIJE(NPT,1)
       IF(XTMP1 >  180.0_SP)THEN
         XTMP = -360.0_SP*TPI+XTMP
       ELSE IF(XTMP1 < -180.0_SP)THEN
         XTMP =  360.0_SP*TPI+XTMP
       END IF	 
!       DLTXE(NPT)=TPI*DLTXE(NPT)*COS(DEG2RAD*YYC)
       DLTXE(NPT)=XTMP*COS(DEG2RAD*Y1_DP)
       DLTYE(NPT)=TPI*DLTYE(NPT)
       CALL ARC(X1_DP,Y1_DP,X2_DP,Y2_DP,DTMP_DP)
       DLTXYE(NPT)=DTMP_DP
       SITAE(NPT)=ATAN2(DLTYE(NPT),DLTXE(NPT))
       DLTXNE(I,1) = DLTXE(NPT)
       DLTYNE(I,1) = DLTYE(NPT)
       DLTXNE(I,2) = 0.0
       DLTYNE(I,2) = 0.0
#      else
       DTMP=DLTXE(NPT)*DLTXE(NPT)+DLTYE(NPT)*DLTYE(NPT)
       DLTXYE(NPT)=SQRT(DTMP)
       SITAE(NPT)=ATAN2(DLTYE(NPT),DLTXE(NPT))
#      endif
     ELSE
       WRITE(IPT,*) 'ISBC(I) NOT CORRECT, I==',I
       CALL PSTOP
     END IF
   ENDDO

   NCV_I  = NCTMP
   NCV    = NCETMP+NCTMP
      
   IF(NCV /= 3*(NT)) THEN
     PRINT*,'NCV IS NOT CORRECT, PLEASE CHECK THE SETUP'
     CALL PSTOP
   END IF
   IF(NCV_I /= 3*N) THEN
     PRINT*,'NCV_I IS NOT CORRECT, PLEASE CHECK THE SETUP'
      CALL PSTOP
   END IF

   DO I=1,NCV_I
     IF(NIEC(I,1) > M .OR. NIEC(I,2) > M)THEN
       write(ipt,*)'problemas',niec(i,1),niec(i,2),m
       CALL PSTOP
     END IF
   END DO

# if defined (SPHERICAL)
   DO I=1, N
     DO J = 1, 3
       J1=J+1-(J+1)/4*3
       J2=J+2-(J+2)/4*3

       DELTX=VX(NV(I,J2))-VX(NV(I,J1))
       DELTY=VY(NV(I,J2))-VY(NV(I,J1))
       X1_DP=VX(NV(I,J1))
       Y1_DP=VY(NV(I,J1))
       X2_DP=VX(NV(I,J2))
       Y2_DP=VY(NV(I,J2))
       CALL ARCX(X1_DP,Y1_DP,X2_DP,Y2_DP,SIDE)
       DELTUX(I,J)= SIDE
       DELTUY(I,J)=TPI*DELTY
       SITAU(I,J)=ATAN2(DELTUY(I,J),DELTUX(I,J))
     ENDDO
   ENDDO

# endif
!==============================================================================!
!  nisbce_1/nisbce_2/nisbce_3:  number of elements with isbce of 1,2,3         !
!  lisbce_1/lisbce_2/lisbce_3:  list of elements with isbce of 1,2,3           !
!  epor                      :  element porosity (=0 if isbce = 2)             !
!==============================================================================!

!
!--COUNT NUMBER OF ELEMENTS OF EACH TYPE (ISBCE=1,2,3)-------------------------!
!
   NISBCE_1 = 0 ; NISBCE_2 = 0 ; NISBCE_3 = 0
   DO I=1,N 
     IF(ISBCE(I) == 1) NISBCE_1 = NISBCE_1 + 1
     IF(ISBCE(I) == 2) NISBCE_2 = NISBCE_2 + 1
     IF(ISBCE(I) == 3) NISBCE_3 = NISBCE_3 + 1
   END DO

!
!--ALLOCATE ELEMENT TYPE ARRAYS LISBCE_1,LISBCE_2,LISBCE_3---------------------!
!
   IF(NISBCE_1 > 0)THEN
     ALLOCATE( LISBCE_1(NISBCE_1) )
   ELSE
!     WRITE(IPT,*)  '!  WARNING               :    NO ELEMENTS WITH ISBCE=1'
   END IF
     
   IF(NISBCE_2 > 0)THEN
     ALLOCATE( LISBCE_2(NISBCE_2) )
   ELSE
!     WRITE(IPT,*)  '!  WARNING               :    NO ELEMENTS WITH ISBCE=2'
   END IF
     
   IF(NISBCE_3 > 0)THEN
     ALLOCATE( LISBCE_3(NISBCE_3) )
   ELSE
!     WRITE(IPT,*)  '!  WARNING               :    NO ELEMENTS WITH ISBCE=3'
   END IF

!
!--LOAD ELEMENT TYPE ARRAYS LISBCE_1,LISBCE_2,LISBCE_3--------------------------!
!
   NISBCE_1 = 0 ; NISBCE_2 = 0 ; NISBCE_3 = 0
   DO I=1,N 
     IF(ISBCE(I) == 1) THEN
       NISBCE_1 = NISBCE_1 + 1
       LISBCE_1(NISBCE_1) = I
     END IF
     IF(ISBCE(I) == 2) THEN
       NISBCE_2 = NISBCE_2 + 1
       LISBCE_2(NISBCE_2) = I
     END IF
     IF(ISBCE(I) == 3) THEN
       NISBCE_3 = NISBCE_3 + 1
       LISBCE_3(NISBCE_3) = I
     END IF
   END DO

!
!--SET ELEMENT POROSITY---------------------------------------------------------!
!
   ALLOCATE(EPOR(0:NT)) ; EPOR = 1.0_SP
   DO I=1,N
     IF(ISBCE(I) == 2) EPOR(I) = 0.0_SP
   END DO
     
   IF (DBG_SET(DBG_LOG))WRITE(IPT,*)  '!  NISBCE/LISBCE/EPOR    :    COMPLETE'
   IF (DBG_SET(DBG_LOG))WRITE(IPT,*)  '!  TRIS/EDGES/CVOLS      :    COMPLETE'

   RETURN
   END SUBROUTINE TRIANGLE_GRID_EDGE
!==============================================================================!
